program chem_thermo_test

   use mod_chem
   use iso_fortran_env, only: dp => real64
   implicit none

   type(type_chem_phases) :: phases
   integer, parameter :: n = 1981
   real(dp), allocatable :: T(:), X(:, :), rho_compared(:)
   real(dp) :: p = 101325.0_dp, rho

   integer :: l, i, ll

   integer :: file_io(100)

   phases = chem_file_load("air_5species_NASA9.yaml", "air")
   allocate (T(n), X(5, n), rho_compared(n))
   open (newunit=file_io(100), file="5S_P_xi.plt")
   read (file_io(100), *)
   read (file_io(100), *)
   do l = 1, n
      read (file_io(100), *) T(l), X(2, l), X(5, l), X(1, l), X(3, l), X(4, l), rho_compared(l)
   end do
   close (file_io(100))

   open (newunit=file_io(1), file="mean_values.plt")
   open (newunit=file_io(2), file="partial_cp.plt")
   open (newunit=file_io(3), file="partial_ethalpies.plt")
   open (newunit=file_io(4), file="partial_entropies.plt")
   open (newunit=file_io(5), file="partial_gibbs.plt")
   open (newunit=file_io(6), file="partial_internal_energies.plt")
   open (newunit=file_io(7), file="partial_cp_mole.plt")
   open (newunit=file_io(8), file="partial_cp_mass.plt")
   open (newunit=file_io(9), file="partial_internal_energies_mole.plt")
   open (newunit=file_io(10), file="partial_internal_energies_mass.plt")

   write (file_io(1), '(1X, A)') "Variables=T, Mw, Temperature, P, rho, delta, cp_mole, cp_mass, h_mass, gamma, soundspeed"

   do ll = 2, 10
      write (file_io(ll), '(1X, A)', advance="no") "Variables=T"
      do i = 1, phases%n_species
         write (file_io(ll), '(1X, A)') phases%species(i)%name
      end do
   end do

   do l = 1, n
      call phases%set_state(T=T(l), X=X(1:5, l), P=p)
      write (file_io(1), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_mean_molecular_weight(), &
         phases%get_temperature(), &
         phases%get_pressure(), &
         phases%get_density(), &
         phases%get_density() - rho_compared(l), &
         phases%get_cp_mole(), &
         phases%get_cp_mass(), &
         phases%get_enthalpy_mass(), &
         phases%get_gamma(), &
         phases%get_sound_speed()
      write (file_io(2), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_cp_R()
      write (file_io(3), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_enthalpies_RT()
      write (file_io(4), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_entropies_R()
      write (file_io(5), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_gibbs_RT()
      write (file_io(6), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_internal_energies_RT()
      write (file_io(7), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_cp_mole()
      write (file_io(8), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_cp_mass()
      write (file_io(9), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_internal_energies_mole()
      write (file_io(10), '(1X, F20.4, 1X, 11G20.5)') T(l), phases%get_partial_internal_energies_mass()
   end do

   do ll = 1, 10
      close (file_io(ll))
   end do

   ! do l = 1, phases%n_species
   !    print*, phases%species(l)%name, phases%species(l)%thermo%cp0_R, phases%species(l)%thermo%dh0_RT
   ! enddo
end
